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ABSTRACT 

We use three-dimensional high-resolution adaptive- mesh-refinement simulations to 
investigate if mechanical feedback from active galactic nucleus jets can halt a massive 
cooling fiow in a galaxy cluster and give rise to a self-regulated accretion cycle. We 
start with a 3 x lO^M© black hole at the centre of a spherical halo with the mass 
of the Virgo cluster. Initially, all the baryons are in a hot intracluster medium in 
hydrostatic equilibrium within the dark matter's gravitational potential. The black 
hole accretes the surrounding gas at the Bondi rate and a fraction of the accretion 
power is returned into the intracluster medium mechanically through the production of 
jets. The accretion, initially slow (~ 2 x lO~^M0yr~^), becomes catastrophic, as the 
gas cools and condenses in the dark matter's potential. Therefore, it cannot prevent 
the cooling catastrophe at the centre of the cluster. However, after this rapid phase, 
where the accretion rate reaches a peak of ~ O.2M0yr~^, the cavities inflated by 
the jets become highly turbulent. The turbulent mixing of the shock- heated gas with 
the rest of the intracluster medium puts a quick end to this short-lived rapid-growth 
phase. After dropping by almost two orders of magnitudes, the black hole accretion 
rate stabilises at ~ 0.006 Mq yr~^, without significant variations for several billions of 
years, indicating that a self-regulated steady-state has been reached. This accretion 
rate corresponds to a negligible increase of the black hole mass over the age of the 
Universe, but is suflicient to create a quasi-equilibrium state in the cluster core. 
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1 INTRODUCTION 



X-ray observations show that active- galactic- nucleus (AGN) 
j ets can inflate large c avities in the hot g a s of galaxy clusters 



(lArnaud et alJ 19841: ICarilli 


et all 11994 


: David et al. 


|200l|: 


Fabian et all l2000l. 20021: 


McNamara et al. l2000l. 


200lh. 


Binnev & Tabor (Il995l) sue 


gested that this injection of me- 



chanical energy could explain why there is little neutral gas 
and star formation at the centre of galaxy clusters, despite 
the fact that the time for radiative co oling is i n ma ny cases 
much shorter than the Hubble time ('Fabian| [T99i ). There 
are three main alternatives to this proposal. The first one 
is that the heating needed to quench a mass ive cooling flow 
is provided by the AGN, but ra diatively (^Ciotti Ostrikerl 
ll997l . l200ll : ISazonov et"al]l2004l ) . and not mechanically. The 
second is that the cluster core is heated from the out - 
side through thermal conduction (e.g. IVoigt Fab "ia^|2004l) . 
The third is that the gas was preheated prior to the epoch 
of cluster formation by e.g. kinetic outbursts at the epoch of 
the formation of elliptica l galaxies and tha t it has slowly 
been cooling since then (iBabul et al.1 l2002l : lOh Bensonl 



l2003l : [McCarthy eral]l2004l ) . However, there are objections 
to both of the flrst two alternatives (see the review by 
"Begelman""2004'), while the work of AGN outflows on the 
intracluster medium (ICM) is observed. 

AGN outflows can be enormously powerful 
(iMcNamara et ^ [2005 1. but the question is how ef- 
fectively this power can be used to heat the ICM in the 
cluster core, e.g. if the kinetic energy of the outflow is not 
converted into thermal energy, a small mass may carry all 
the energy and escape with it from the cluster. Moreover, 
luminous AGNs are short-lived. It is not clear that such 
erratic energy sources will have any long lasting eflect. 
The problem is complicated because AGN outflows are 
non- spherical, time-dependent and turbulent, so the only 
proper method of modelling the interaction between the 
jets and the ICM is numerically. 



In the last few years, a number of groups have 
made an intense effort to simulate the impact of 
radio jets and bubb le s on the s t ructu r e of the ICM 



(IChurazov et al. 2001: iQuilis et al. 


I2OOII: iRevnolds et al. 


I2OOIL I2OO2I: iBasson k Alexander 


I2OO3I: lOmma et al. 
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20041: [Omma Binneyl l2004l: iRuszkowsk i et alJ l2004l : 
Siiacki Springell 1200^1 Vernaleo k Reynol ds 2006). This 
work has given us a reasonable picture of how the lobes 
of radio sources are formed and of how they evolve after 
the central engine has stopped pumping energy into them. 
It has helped clarifying the mechanisms through which 
the mechanical energy can be thermalized, and the issues 
involved, e.g. the viscosity of the ICM, and the difference 
between heavy and light jets. 

In all this studies, the power and the duration of the 
AGN phase were set by hand. The AGN was able to affect 
the ICM, but the change in the ICM properties was not al- 
lowed to have any feedback on the AGN power. Here, we 
remove this limitation by considering a model in which: i) 
the jet power is proportional to the accretion rate of the 
central black hole, and ii) the black hole accretion rate de- 
pends on the central density and temperature of the ICM 
through the iBondil (| 19521 ) model. Thus, we close the loop 
and investigate the self- regulation of the AGN - cooling flow 
system directly. 

In this paper we start with a 3 x 10^ Mq black hole at 
the centre of a 1.5 x 10^^ M© halo. We simulate how this 
system evolves due to radiative cooling and AGN feedback 
in order to establish the pattern of black hole accretion and 
the capacity of the system to approach a state of equilibrium. 

The paper is structured as follows. In Section 2 we de- 
scribe our model for the hydrodynamics of the ICM and 
its interaction with the central black hole. In Section 3 we 
present the simulations and their results. In Section 4 we 
discuss the conclusions that can be derived from these re- 
sults. 



2 THE MODEL 

2.1 The intracluster medium 

We assume that we can treat the ICM as an ideal gas and we 
follow its dynamics in the gravitational potential of a static 
spherical dark matter halo, the radial density pro file of whic h 
is the described by the NFW model ( Navarro et al.lll997l l. 
We include energy dissipation due to radiative cooling as 
well as the injection of mass, momentum and energy due to 
the presence of an accreting black hole at the centre of the 
gravitational potential. The equations of motions in their 
conservative form are: 



d 

— (pv) + V • (pV V) + Vp = -pV(/)dm + '^Pjll^ 

ot 



(1) 



(2) 



^(pe) + V. 



pv 
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-PV • Vflidm - C + IVIPj (3) 



(4) 



where p is the density, v is the velocity, e is the specific 
total energy, p is the pressure, (^dm is the gravitational po- 
tential of the dark matter, C is the radiated power per unit 
volume), 7 is the ratio between the specific heats at con- 
stant pressure and volume, Mjet, Pjet and Pjet are the mass, 
momentum and energy injection rates due to feedback from 
black hole accretion, is a function that specifies the spatial 



distribution of the injection (the integral of IV^I over space 
is equal to unity) , while is the unit vector in the positive 
direction of z axis, chosen to coincide with the axis of the 
jets coming out of the central source. 

The radiated power per unit volume has the form 
C = nHA(T, Z), where nu — fup/rup is the number den- 
sity of hydrogen nuclei. Here rup is the proton's mass and 
/h = 0.75 is the mass fraction of the ICM in hydrogen. The 
cooUng function A only depends on the temperature T and 
the metallicity Z. It is computed fro m the collisional equi- 
librium cooling functions tabulated in [Sutherland Dopital 
(|l993l). The temperature from which we interpolate the cool- 
ing rate is computed from p and p using the equation of 
state for an ideal gas with an atomic weight of /i :^ 0.6mp, 
the value for a totally ionised plasma with 25% helium in 
mass. All our calculations are for a metal abundance of one 
third the Solar value. 

In our simulations, we assume that all the baryons are 
in the ICM and that the cosmic baryonic fraction is 10%. 
We neglect the fact that the ICM has a non-zero total an- 
gular momentum. We assume that initially the ICM is in 
hydrostatic equilibrium within a static NFW potential. Ap- 
pendix A contains the technical details of how we set up our 
hydrostatic initial conditions. 



2.2 Black hole accretion and feedback 

We start with a supermassive black hole at the centre of 
the gravitational potential. The bl ack ho l e grow th rate is 
computed with the Bondi formula. iBondil (|l952[ ) studied a 
spherical stationary solution where the gas is static and ho- 
mogeneous at infinity and in free fall close to the black hole. 
The transition between the two regimes occurs at the sonic 
radius rs, where the circular velocity on a Keplerian orbit 
around the black hole is equal to the sound speed Cg of the 
ICM (GM./rs = Cs). Under this hypothesis, the black hole 
accretion rate is 



M« — 47rrs pCs 



:47r(GM.)^4 



(5) 



where p is the gas density at the sonic radius. In reality, the 
gas distribution around the black hole is non-spherical and 
turbulent, so it is difficult to define the sonic radius. For this 
reason, we adopt a practical approach and determine p and 
Cs as the density and the volume-weighted sound speed in a 
fixed small but resolved volume around the black hole. It is 



also worth noting that, for 7 = 5/3, pc^ 



(5/3) 



3/2-3/2 

where s = pp~^ — sq exp[(7 — l)pmpk~^ (S — So)] is a mea- 
sure of the entropy per unit mass, of the ICM {k is the 
Boltzmann constant and so is the value of s that corresponds 
to the zero point of the specific entropy, 6*0) • Hence, in the 
Bondi model, M« depends on the state of the ICM through 
its entropy only. 

We assume that black hole accretion is accompanied 
by the acceleration of jets. Jets emit radio synchrotron 
radiation because they are magnetised. A magnetic field 
of 1 nT in vacuum corresponds to a magnetic pressure of 
^ 10~^^dyncm~^, comparable to the thermal pressure of 
the ICM. Magnetic fields should therefore be important. 
They are also believed to play a role in collimating the jet 
itself. However, we assume here that we can model jets using 
only Euler's equations. This approximation is valid if one is 
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interested in the dynamics of the cocoons inflated by the 
jets, rather the dynamics of the jets themselves. These hot 
bubbles are responsible for the mechanical work on the ICM. 
As long as the energy and the momentum pumped into the 
jets are calibrated correctly, one will end up reproducing the 
expansion of the cocoons appropriately. 

We assume that the power of the relativistic particles 
pumped into the jets is proportional to the accretion rate of 
the black hole and that on subparsec scales, where jets are 
accelerated, the outflow rate is equal to the accretion rate 
and all the energy is kinetic. With these hypotheses, the 
rates at which momentum and energy are deposited into 
the jets are: 

Pj V2^ M.c (6) 

Pj = e M.c^ (7) 

Here c is the speed of light and e is the efficiency with which 
the rest-mass energy of matter accreted onto the black hole 
is used to power the jets. 

As jets travel from the subresolution region where they 
are accelerated to the smallest scale that we can resolve 
in our simulations, they entrain additional material from 
the interstellar medium of the black hole's host galaxy. This 
increases the mass outflow rate, while conserving both the 
momentum and the total energy. Defining 77 as the jet mass 
loading factor, we have: 

Mj = ry M. (8) 

This process occurs at scales far below our resolution limit. 
We consider here only the large scale jets, for which the 
velocity is much smaller than the speed of light. The jets are 
launched in two opposite directions from the central ^ 3 kpc 
at a speed — p-J M-^ — 77"^ Y^2e)c. Note that this speed is 
consistent with momentum conservation in the jet. A mass 
loading factor of 77 = 1 corresponds to the case without 
entrainment, in which all the energy is kinetic. For 77 > 1, 
Pj > l/2Mjt'j^. The difference between the two is added as 
thermal energy to the gas injected at the base of the jets. 
The physical reason is that the entrained gas is shocked. 
Therefore, kinetic energy is converted into heat. 

In our mathematical formulation, Mj, and Pj are 
function of time only and as such they do not contain any 
spatial information, which is, instead, in the function '0 
(Eqs. [T][2l[3|. The precise form of the function -0 is dictated 
by computational reasons and has no physical origin. Its con- 
sideration is therefore postponed to the section that deals 
with the numerical aspects. 

2.3 Choice of physical parameters 

Our choice of environmental and black hole parameters is 
dictated by observations of M87, although our study is not 
an attempt at realistically modelling any individual object. 
In particular, we have neglected the stellar potential, which 
dominates the gravity in the central parts of M87. The halo 
has a mass of Mvir = 1.5 X I^^^Mq. The virial radius, 
Tvir — 1.4 Mpc, is computed by assuming that the virial 
density is Ac = 101 times the critical cosmic density. The 
core radius of the dark matter distribution, ro — 250 kpc, 
is computed from the virial radius with the concentration 
parameter given by Eq. ()A2p . 



Table 1. Simulation parameters (halo, black hole, numerical) 



Mvir 


1 c 1 n 14 A /r 
i.5 X iU Mq 


^vir 


1.4 Mpc 


^0 


250 kpc 


M. 


3 X lO^M© 


e 


0.1 


V 


100 




3.2 kpc 


h 


2.5 kpc 


^box 


648 kpc 


Ar 


0.64 kpc 


^sim 


12Gyr 



The initial mass of the central black hole is M« = 3 x 
10^ M0. It is commonly assumed that the accretion of matter 
onto a black hole releases energy at ^ 10% efficiency, so that 
the accretion power is Pacer ^ 0.1M«c^, and that most of 
this power is emitted as light. This is almost certainly true 
for the AGNs detected in optical surveys. In fact, there is 
reasonably good agreement between the cosmic black hole 
density inferred from mass estimates in the local Universe 
and the density derived by integrating the optical and X- 
ray emission from AGrNs over the entire l ife of the Universe 

(|Yu Tremainell2002l : lBarger et al.ll2005l ). 

However, both observation (M87; IPi Matteo et al.l 

2OO3II and theory (ADIOS ad iabatic inflow-outflow solution; 
Blandford BegelmanI [TqQqI ) suggest that this may not be 
true when <C MEdd, where MEdd is the accretion rate 
needed for the AGN to radiate at the Eddington luminosity. 
In such cases, the most probable outcome is that most of 
the power is released mechanically. 

In the simulation that we present in this article, the 
maximum accretion rate is M« ^ O.2M0yr~^ (Section 3). 
For M. = 3 X 10^ Mo, this gives M. - 3 x 10"^MEdd. 
Therefore, we are at all times inside the regime described by 
ADIOS, where we can assume that all the power generated 
by the AGN is channelled into the jets. If matter accreted 
by the black hole releases energy at the canonical efficiency 
of ~ 10%, this gives us a canonical value of e ^ 0.1. 

We have set the mass loading factor to 77 = 100 so 
that plasma is injected away from them black hole at v-^ ^ 
1350 km s~^ in agreement with observations of Fanaroff- 
Riley I sources. This speed is much lower than the one at 
which particles are accelerated close to the central engine be- 
cause the jets have entrained mass and have therefore slowed 
down before reaching the 3 kpc scale at which we introduce 
them in our simulations. The parameter 77 also determines 
the temperature of the plasma, since the rates at which ki- 
netic energy and thermal energy are produced must add to 
the total power that the AGN deposits into the IGM. Ghoos- 
ing ^ 1350 km s~^ implies that the kinetic energy is only a 
small fraction of the total energy of the jets . Our jet model 
is the refore intermediate between those of iRevnolds et al.l 
(|2002l ) and lOmma et all (|2004l l although closer to the for- 
mer. This non-relativistic speed also ensures that we are 
allowed to model the system with classical hydrodynamics. 
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2.4 Numerical aspects 

The adapt ive mesh refin ement (AMR) hydrodynamic code 
RAMSES (|Tevssiei]l206i ) is used to integrate the equations 
of the conservation of mass, momentum and energy on a 
three-dimensional Cartesian grid. The size of a computa- 
tional cell is 2~^Lbox, where Lbox is the size of the cubic 
computational box and I is the level of refinement. In AMR 
codes such as RAMSES, I can have a large value where high 
spatial resolution is needed (e.g. in the cluster core or along 
the jet axis) without slowing down the computation by im- 
posing high resolution over the entire box. The time step 
is also adaptive and is determined for each level of refine- 
ment independently by using standard stability constraints 
for hydrodynamic solvers. 

Our refinement strategy is based on the geometry of the 
problem, which is a combination of the spherical geometry 
of the potential and the cyUndrical geometry of the jets. 
The computational domain is a cube of size Lbox — 650 kpc, 
subdivided in large cells of size 2~^"^^"Lbox, where ^min = 5. 
We define an ellipsoidal region, centred on the cluster core, 
with minor axis a ^ 40 kpc and major axis b ^ 160 kpc. 
The major axis is aligned with the jet axis. This region is 
refined up to the maximum level, imax = 10. The cell size 
at the highest level of refinement is Ar ^ 0.6kpc. We then 
progressively de-refine the grid outside this central region 
in order to save memory and computing time. We obtain 
a total number of AMR cells Ntot — 5 x 10^, so that one 
run takes approximately 100 hours wall clock time on 32 
processors for a total number of time steps roughly equal to 
10^ 

Jet simulations are very time consuming because the 
time step, controlled by the Courant stability condition, is 
very small. This has prevented us from performing higher 
resolution simulations to estimate the convergence proper- 
ties of our settings. In cooling- flow simulations, the gas tem- 
perature is lower and the number of time steps is one order of 
magnitude smaller. We have therefore been able to increase 
the size of the box up to twice the virial radius and to raise 
the number of cells up to 5 x 10^. In this case, we obtain a 
cooled gas fraction of ^ 7% instead of the one obtained in 
our 'fiducial' run (the one with the same resolution and box 
size that we use for jet simulations), which is around 3.5% of 
the total gas mass. Due to a limited resolution and box size, 
we might therefore underestimate cooling by a factor of ~ 2. 
However, due to the self-regulating nature of AGN feedback 
(Section 3), we are confident that our ultimate conclusions 
remain unaffected by these numerical issues. 

The central region that we use to compute the accretion 
rate of the black hole is a small cylinder aligned with the jet 
axis. Its radius and half-height are rj = 3.2 kpc and /ij = 

2.5 kpc, corresponding to 5Ar and 4Ar, respectively. The 
same cylinder is used to define the spatial distribution of 
the jet injection of mass, momentum and energy, according 
to 

1 / x^±y^\ z 

^=2^^"P^-^^JXf 

for + 2/^ < and \z\ < /ij, and ip — everywhere 
else. This scheme for launching jets was introduced by 
lOmma et al.l (|2004 ). We have verified that, for the same 
setup, our results and theirs coincide. 



The length-scale of our cylinder is much larger than 
the estimat e of ^ 120 pc for the Bondi radius in M87 by 
[ah en et al. 1 (|2006"). One could therefore reasonably worry 
that smoothing over such a large region may lead to an un- 
derestimate of the actual mass inflow rate at the Bondi ra- 
dius. However, since the nearly flat density and temperature 
profiles that we measure in the central part of our simulated 
cluster (see Figure [2]) indicate that the accretion rate com- 
puted with Eq. (|5]) should not be highly sensitive to the 
radius at which the density and temperature are evaluated. 



3 SIMULATIONS AND RESULTS 

In the course of this study, we have run three simulations. 
The first one is a purely hydrostatic simulation, where we 
have switched off cooling and feedback: we have verified that 
the density distribution of the ICM remains constant for the 
whole duration of the run, up to tsim 12 Gyr. As this first 
simulation is just a check on the accuracy of our hydro- 
static equilibrium, we shall only comment on the other two. 
The second simulation is a standard cooling- flow simulation, 
where we have included cooling, but without AGN feedback. 
The third simulation is our main result: both cooling and 
AGN feedback are considered. The global model parameters 
are summarised in Table 1 and are the same for all three 
simulations. 

In the standard cooling-flow simulation, only at t ~ 

4 Gyr does the gas start to cool significantly. This is due to 
the initially high entropy of the core of the simulated clus- 
ter (e.g. Oh & Ben son! l2003l : iMcCarthv et al.H2004 ). How- 
ever, as it starts condensing to the centre, cooling becomes 
catastrophic (in the bremsstrahlung regime, the cooling rate 
varies with the square root of the temperature, but the 
square of the density) . At t 6 Gyr a dense cool core is 
clearly visible (Fig. la). The temperature in the central 6 kpc 
has decreased by a factor of ^ 4 and the density has gone 
up by a factor of ^ 10 (Fig. 2). The entropy has dimin- 
ished accordingly. In Fig. 2, we have shown the entropy 
computed as Tnjj~^, where T is the temperature in keV 
and riH is the number density of hydrogen nuclei in cm~^. 
This definition is related to the one given in Section 2.2 
by s = kT{nnlJim^Y~^ . By the end of the simulation, at 
t = 12 Gyr, the central density is more than a hundred times 
higher than its initial value, the central temperature and 
entropy have decreased dramatically, while the pressure has 
grown by a factor of ^ 2. The pressure has changed more 
moderately than the other three quantities because the cool- 
ing time is longer than the dynamical time and the ICM has 
cooled quasi-statically. 

The gas in our simulations is either hot, with a tem- 
perature of the order of the virial temperature of the halo 
(T - 10'^ K), or cold (T - 10^ K), with a negligible fraction 
of the mass at intermediate values. It is therefore straight- 
forward to separate the cold gas from the hot gas and to 
see that a large mass of cold gas builds up in the central 
few kiloparsecs. This mass is of - 10^° Mo at t - 6 Gyr. 
By the end of the simulation, it grows above 4 x 10^^ M© 
(Fig. 3). This cannot be appreciated from Fig. la because 
gas that cools below 10^ K goes down to 10^ K very rapidly, 
fragments into clumps and it gets squeezed into a few dozens 
of cells. They contain all the cold gas and occupy a central 
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b. Model with AGN feedback 
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Figure 1. Maps of the emission- weighted temperature along the hne of sight in the coohng flow simulation (flrst row) and in the feedback 
simulation (second row and third row). The initial condition (t = 0) is the same for both simulations. Each map corresponds to a square 
fleld with a side of 260 kpc. The axis of the jets is perpendicular to the line of sight and coincides with the vertical direction on the maps. 
The temperature scale used to colour the snapshots is logarithmic. 



region with a radius of ^ 4 kpc. The lack of net angular 
momentum in our initial gas distribution, the absence of a 
prescription for star formation and supernova feedback and 
the limited number of cells in this central region imply that 
we have evolved our cluster to the limit where our physical 
description breaks down. 

A resolution study shows that the reason why only 
~ 3.5% of the baryons cool is that our simulation is not 
fully converged but achieving the convergence resolution is 
impossible with the short timesteps of the jet simulation. 
We have therefore shown the outputs of the cooling flow 
simulation at the same resolution used for the feedback sim- 



ulation. We shall discuss later how this is likely to affect our 
conclusions. 

In the simulation with feedback, the interaction of the 
black hole with the ICM starts gently because the initial 
accretion rate, M, ^ 5 x lO^^M© yr~\ is low (Fig. 4). This 
is due to the high entropy of the cluster core in the initial 
condition. The AGN empties a narrow tube along the jet 
axis (Fig. lb and Fig. 5). The jets look like spring onions, 
with a thin round part at the bottom, which continues into a 
stem. This can be explained by the pressure of the jets when 
they come out of the injection zone, which is slightly lower 
than the pressure in the cluster core. So the jets are squeezed 
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Figure 2. Density, temperature, pressure and entropy profiles in the cooling fiow model (red) and in the model with AGN feedback 
(black) at t ~ Gyr (solid lines), t ~ 6Gyr (dotted lines), t ^ 8 Gyr (dashed lines, only shown for the model with AGN feedback) 
and t 12 Gyr (dashed-dotted lines). The initial condition is the same for both the cooling fiow and the AGN feedback simulation. 
The density nn is the number of hydrogen atoms divided by volume. T is the emission-weighted temperature. The 'entropy' Tn^^ is 
weighted over the X-ray emission, too The pressure p is the volume-weighted average. All four quantities are calculated in concentrical 
spherical shells and are given as function of the radius r of a shell from the centre of the cluster. 



at the sides, until they reach the pressure of the ICM. Due to 
the continuity equation, the shrinking of the pipe's section 
is accompanied by an increase of the flow speed. This is seen 
as an increase of the ram pressure pv^ at the points where 
the jets on the temperature maps become narrower (Fig. 5). 

Using directly the density and temperature maps, we 
can determine that the jets are propagating into the cluster 
core at a subsonic speed of about 500kms~^. This subsonic 
propagation can also be demonstrated through the absence 
of hotspots in the lobes. The temperature of the plasma in 
the jets has a maximum at the central source and decreases 
moving outwards. This is consistent with the notion that 
we are simulating a Fanar off- Riley type I or core-dominated 
radio source. The behaviour of the density mirrors that of 
the temperature. The density grows outwards along the jets, 
as the outflowing plasma slows down and gathers in the jet 
lobes, where it mixes with entrained material. 

At early times, the AGN has little effect on the state 
and evolution of the ICM in the cluster core. The gas en- 
tropy in the cluster core is rather high, resulting in a rather 
low black hole accretion rate of M, < 10~^Mq yr~^. There- 



fore, even with the high efficiency we have considered in the 
accretion-ejection conversion process, the power injected by 
the AGN into the ICM is less than the rate at which the 
ICM looses energy through radiative cooling (Fig. 4). Still, 
the injected energy delays the onset of the fast cooling phase 
from t ^ 4 Gyr up to t ~ 6 Gyr (see Figs. 3). Another reason 
for this low overall efficiency at early time is that the high 
entropy plasma in the jets flows outwards, rather than en- 
riching the entropy of the ICM around the black hole. This 
can be seen by comparing the cooling flow simulation and 
the AGN feedback simulation using Figure lat t = 6 Gyr. 
Both contain a cool core at T < 10^ K, but the most vis- 
ible difference is at r > 100 kpc, where the ICM has been 
heated more efficiently by the AGN feedback than in the 
pure coofing flow simulation. 

By t ^ 6.5 Gyr, the cooling catastrophe, although it 
was delayed by almost 2 Gyr, still occurs. As the gas in 
the cluster core is cooling and condensing faster and faster, 
the accretion rate of the black hole grows by two orders 
of magnitudes in ^ 0.5 Gyr, and reaches a peak of M, ^ 
O.IM© yr"^ at t 7 Gyr. At this rate, the jet power, Pj™^"" - 
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Figure 3. The mass of the cold gas (T < 10^ K) in the computa- 
tional volume as a function of time. The dotted line corresponds 
to the cooling flow simulation and the solid line to the feedback 
simulation. 



Figure 4. The accretion rate of the black hole, M», and the cor- 
responding jet mechanical luminosity, 0.1M«c^ in the simulation 
with AGN feedback (solid line). The dotted line and the dashed 
line show the X-ray luminosity emitted by the ICM in the cooling 
flow simulation and the feedback simulation, respectively. 



lO^^ergs"^, exceeds the bremsstrahlung luminosity of the 
cluster by a factor of ^ 10^. The thermal pressure p at the 
base of the jets and the ram pressure pv^ throughout the 
length of the jets are no longer negligible compared to the 
pressure of the ICM in which the jets propagate. 

Given the new physical state in the cluster core, the 
jet propagation is no longer stable. The Kelvin-Helmoltz in- 
stability sets in and rapidly amplify small scale fluctuations 
already present in the flow. This effect is so strong that 
the flow at t = 6 Gyr is no longer axisymmetric, although 
our initial conditions were perfectly axisymmetric (Fig. 5). 
A complete analytical description of jet instabilities is be- 
yond the scope of this paper fsee Smith et al. (1983) for a 
comprehensive study). Our numerical experiment suggests 
that these instabilities grow faster in the phase of rapid ac- 
cretion: the jets become turbulent and break up into bub- 
bles, which, because of their lower density than the ICM, are 
buoyantly driven outside the cluster core. In this turbulent 
atmosphere, bubbles may rise in directions not aligned with 
the jets' axis. Fig. 6 shows such an extreme situation: the jet 
pointing downwards blows a bubble that rises perpendicular 
to the jet. The hydrodynamic interaction of the jets with the 
ICM is, therefore, a complicated process with many different 
aspects to it. The jets inflate cavities and through their ex- 
pansion blow the gas in the centre of the cluster away. They 
heat gas through shocks, but also through turbulent mixing 
of jet material with the ICM. AGN feedback reheats most of 
the cold gas that has accumulated at the centre of the clus- 
ter during the cooling catastrophe and causes M, to drop 
by a factor of '-^ 50 at 7 Gyr < t < 8 Gyr. After reaching its 
maximum value of M« ~ 2 x lO~^M0yr~^, the accretion 
curve is reasonably modelled by a decreasing exponential, 
with an e-folding time scale of tufe ^ 2.5 x 10^ yr, consistent 
with the typical life time of a radio source. 

At t > 8 Gyr, the central temperature of the hot gas, the 
mass of the cold gas and the black hole's accretion rate settle 
respectively to T - 3 x lO'^ K, Mcoid - 5 x lO'^M© and M. - 
2 X 10~^, without substantial variations until the end of the 
simulation at t 12 Gyr (Figs. 2, 3 and 4, respectively). 
This strongly suggests that an equilibrium state has been 



reached. This new equilibrium regime can be characterised 
by a highly turbulent and highly pressurised state within the 
cluster core. These new properties play an important role 
in the system's self-regulation. Before the cooling flow goes 
unstable, the hot plasma flowing out of the central source 
easily escapes the cluster core, funnelled through the high 
density ICM by a straight and empty channel excavated by 
Kelvin-Helmholtz stable jets. In the new equilibrium state, 
this is no longer the case: th e jets must work constantly 
to open a path for themselves. lOmma BinnevI (|2004l l had 
also pointed out the importance of closing the channel for 
effective feedback. Moreover, the higher pressure of the ICM 
contributes to confine the outflow within the central region 
and force most of the mechanical energy to remain inside the 
core. Therefore, the jets at t = 12 Gyr are slower, hotter and 
ultimately more effective than those at t = 6 Gyr (Fig. 5): a 
higher fraction of the AGN power is converted into heat that 
compensate bremsstrahlung cooling in the central region. 
This explain why the inner region of the cluster has reached 
this new hydrostatic and quasi-adiabatic regime. 

In the catastrophic cooling phase (t ^ 6 Gyr), the en- 
tropy of the ICM within the core decreases from 100 keV cm^ 
to 20 keV cm^. At t 8 Gyr, after the strong AGN outburst, 
the central entropy is nearly back to its initial value. By 
t ~ 12 Gyr, the entropy of the cluster core has started to 
decrease again, albeit by a small amount and over a rather 
long period of several Gyr. This slow evolution is due to a 
secular rise of the central density, as the central tempera- 
ture remains quasi constant (Fig. 2). In Section 4 we shall 
argue that the temperature does not change much because 
it adjusts itself to the value at which AGN heating balances 
radiative cooling in the cluster core. We should make note of 
the fact that the profiles in Fig. 2 are useful to understand 
the physical processes in the simulated cluster but cannot 
be taken as accurate predictions to be compared with ob- 
servations because of the idealised gravitational potential 
used in the simulations, which neglects the dominant stel- 
lar component at < lOkpc fe.g. iMamon Lokasll2005l l. We 
have also neglected the role of cosmic rays as an additional 
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Figure 5. The thermal pressure, p = pkT / iirrip^ the ratio of the ram pressure to the thermal pressure, pv^ jp^ and the temperature, T, 
of the gas at two different times = 6 Gyr in the upper row and t = 12 Gyr in the lower one) on a section that contains the jet axis and 
is orthogonal to the line of sight used to plot Fig. 1. The charts have a side of 260 kpc. All colour scales are logarithmic. 



non-equilibrium component, whose observational signature 
mainly resides in the radio emission of galaxy clusters. 

The cooling of gas in the outer parts of the cluster core 
causes the cluster gas to flow inwards and slowly compress 
the central region (Fig. 1). This secular compression is con- 
trolled by the quasi-adiabatic equilibrium we have reached 
within the cluster core. In a more realistic cosmological envi- 
ronment, the merging of clumpy satellites of cold gas would 
have triggered another phase of strong AGN activity, fol- 
lowed by another long period of quasi-equilibrium. 



4 DISCUSSION AND CONCLUSION 

Our views on the thermal evolution of the hot gas in 
galaxy clusters have changed considerably in recent years 
due to data from CHANDRA and XMM ("Fabia n et all 
These X-ray satellites have found no evid ence for 
gas below one third of the virial temperature (| Allen et al.l 
I2OOII : IPeterson et al.|[2003h . Millimetre observations of CO 
lines with the IRAM 30 Metre Telescope, the James Clerk 
Maxwell Telescope and the Caltech Submillimeter Ob- 
servatory have det ected the elusive cold molecular gas 
(|Salome et al.ll2006l and references therein). However, mea- 
sured cooling rates are sub stantially lower than in ferred from 
X-ray luminosity (also see iBregman et al. I l2006h . Inefficient 
cooling in massive haloes is also necessary to explain several 



galaxy properties ( Cattaneo et al. I I2OO6I . I2OO7I : ICroton" et al. 
I2OO6I : iBower et allbood l' 



There is increasing consensus that AGN feedback is 
the reason why th e cooling rate is much lower than one 
would expect (e.g. iTabor Binney _[l993': Bin nev h Taborl 



Reynolds et al.l 
Begelman 

2004: IR uszkows ki et al 



2004 



1995; iTucker DavidI Il997l: "Tciotti Ostrikerl 
I2OOII : IChurazov etHI I2OOII : lOuilis et al ' 



2OQIL I2OO2I: Sasson h Alexand er! |200^ 
lOmma et al. 2004; Omma Binnevl 




Brighenti Mathew; 



Mathews et al II2OO6I: iN^ sser et J] l2006l : ISiiacki h Sprint 



I2OO4] : 

bool 



IZanni et al.l 
iFabian et al.l 



2005 



2006 



20061 : [ Vernaleo &: Revnolds» i2006t ) . Non- gravitational heat 



ing is also needed to explain the high entropy floor of 
the ICM and why the observed relation between X-ray 
temperature a nd X-ray luminosity d e viates from the virial 
relation (e.g. McCarthy et al.l I2OO2I : iRovchowdhurv et al.l 
I2OO4I : iBorgani et al.ll2005l l. AGN feedback can be radiative 
and mechanic. Most of the computer simulations that 
follow the hydrodynamics of the interaction of AGNs with 
the ICM have concentrated on mechanic feedback because 
of the direct observational evidenc e that jets from AGNs 
can open large cavities in the ICM (|McNamara et al.]|2005l : 
Fabian et al. 2006). 

However, while substantial work has been put into un- 
derstanding the mechanisms through which jets affect the 
ICM, the exploration of the long-term outcome of this in- 
teraction and of the feedback that it returns to the central 
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engine has just started. lOmma Binne 3 (|2004h compared 
two simulations of a cooling flow cluster in which jets were 
turned on at different times. They assumed that a later ac- 
tivation of the central source results in more powerful jets, 
since the central density of the ICM is higher at a more 
advanced stage of the cooling catastrophe. They used this 
as an argument to suggest that the interaction of the black 
hole with the ICM is structurally stable. However, in nearly 
all the studies done until now, the power and the duration 
of the jets were set by hand following observational esti- 



mates. 



Brigh enti Math ewsl (|2006[ )'s 2D simulations and 



IVernaleo Revnolds (200^) 's 3D study, whose simulations 
did not last long enough to attain the self-regulated regime 
described in this article, are the only ones so far who have 
attempted to relate the properties of jets to the cooling rate 
of the ICM. 

This paper fills this gap in our comprehension of the 
problem by closing the feedback loop with a self- consistent 
model of the interaction between the black hole and the 
ICM. The accretion rate of the black hole is computed from 
the density and temperature of the ICM. The power of the 
jets is proportional to the accretion rate of the black hole 
and the jets affect back the accretion rate of the black hole 
by changing the properties of the ICM. 

The specific model used to calculate the acc r etion rate 
of the black hole is the one proposed by Bondi (1952). In 
this model, which corresponds to spherical accretion of gas 
that is supported by thermal pressure and static at infinity, 
the accretion rate of the black hole is M» oc M^s~^^'^ where 
s oc pp~^^^ oc Tn^^^ is a measure of the entropy of the gas 
around the black hole (Section 2.2). 

In this work, we have simulated the evolution of the 
ICM in a halo of Mh aio 1-5 x ^^I'^Mc, described by the 
NFW density profile ([Navarro et al.lll99 ?) and containing a 
central black hole of M, = 3 x lO^M© starting from hydro- 
static initial conditions. In our simulations, the growth of 
the black hole is so small that M« is nearly constant, but 
its value is important because it determines the power of 
the central source. The black hole feeds energy back to the 
ICM by the injection of hot plasma with an outfiow speed 
of V] ^ 1400 km s~^) in two symmetric volumes at a small 
distance from the black hole. 

From our simulations of the interaction of the black 
hole with the ICM, we learn two new fundamental results. 
In presence of a massive black hole, a cooling-fiow cluster 
switches on an AGN of ^ 10^^ ergs~"'^soon after the onset of 
the cooling catastrophe. Its lifetime is short, 2.5 x 10^ yr, 
because twin jets blow away and reheat the gas in the cluster 
core. After the AGN phase, there are at least 4 Gyr during 
which the accretion rate of the black hole and the mass of 
the cold gas in the cluster core stay constant, suggesting 
that a self- regulated equilibrium has been reached. 

Turbulence stirred by jet instabilities mixes the jet 
plasma w ith the ICM in the c lu ster core, a s previ ously 
found bv IChurazov et aD (|200lh . lOuilis et all (|200lh and 
iRevnolds et al.l (|2002h . A highlv turbulent and highly pres- 
surised state of the ICM is essential for the effectiveness of 
our feedback mechanism. Without it, the jet energy would 
have left the cluster core and the AGN would have not 
heated the cooling gas efficiently. Before the jets grew un- 
stable, most of the power was indeed deposited far from the 
centre. Self-regulation was then impossible. 



Figure 6. The cavities that AGN feedback has opened in the 
hot gas appear as lower density regions. Above is a density map 
extracted from the simulation at f ~ 8 Gyr. Below is an observa- 
tional image of the radio emission from the synchrotron plasma 
that fills the X-ray cavities in M87 (fOwen et al.l l2000V The two 
figures have the same physical size (~ 70kpc x 70kpc) and show 
an interesting morphological similarity. 



The violent outburst following the cooling catastrophe 
causes the thermal and turbulent pressure to rise signifi- 
cantly, until most of the jet energy is confined inside the 
cluster core. This confinement allows energy released by the 
black hole to have an immediate feedback on its accretion 
rate. This ensures that an equilibrium is reached between 
the accretion rate of the black hole and the temperature of 
the surrounding gas. 

The equilibrium state can be determined by balanc- 
ing AGN heating with radiative cooling in the cluster core, 
assuming that it is now a closed system, so that all the 
energy released by the AGN remains within the bound- 
aries of the core. Eqs. (5) and (8) give the heating rate 
Pj '2^ 47rec^(GM,)^mpnHc7^ oc nnT"^/^, while the total ra- 
diated power is given by Lx = nHA(T)MH/mp oc riYiT^^^ ^ 
where Mh, the total hydrogen mass in the core, is assumed 
constant (Mr ~ 10^^ M©). Given Mr and M«, there is one 
equilibrium temperature for which Pj = Lx and its value its 
independent of nn- From this simple argument, we compute 



10 A. Cattaneo, R. Teyssier 



T = 2 — SxlO^K. This figure agrees with the results of our 
simulations (Fig. 2) and is slightly higher than the virial 
temperature of our simulated cluster, Tvir — 1.7 x 10^ K. 

If the mass of the black hole had been significantly lower 
than the value considered here, the equilibrium temperature 
would have been much lower than the virial temperature. As, 
by definition, the virial temperature is that at which ther- 
mal pressure balances the gravitational force, AGN heating 
would not have been able to stop the contraction of the clus- 
ter core with such black hole mass. The black hole would 
have had to grow considerably for feedback eventually to 
become efficient. 

In addition, if the local cooling rate exceeds the local 
heating rate around the black hole, the net effect is cool- 
ing, the temperature decreases, the density increases and 
the heating rate goes up. If, on the other hand, the heating 
rate becomes larger than the cooling rate, than the temper- 
ature increases, the density decreases and the heating rate 
goes down. This thermal equilibrium is, therefore, likely to 
remain stable and self-regulated for several Gyr, as our sim- 
ulations have demonstrated. The infall of cold gas clumps 
from outside the c entral region (such as those observed by 
ISalome et al.ll2006l ) may eventually give rise to a new period 
of AGN activity. 

The simulations presented in this paper prove that the 
interaction of the AGN with the cooling flow evolves natu- 
rally towards a self-regulated equilibrium state after a short 
phase of strong activity. We suspect that this feature is ro- 
bust and does not strongly depends on the mechanism that 
thermalizes the power channelled into the jets. 

Our simulations do not provide an accurate description 
of an individual object such as M87, since, among other 
approximations, we have neglected the stellar contribution, 
which dominates the gravity in the central few kiloparsecs. 
However, the halo and the black hole parameters were chosen 
to mimic the values appropriate for M87. Therefore, a careful 
comparison with M87 is useful because it gives us an idea 
to what extent this simple model represents the physical 
reality. 

iBoh ringer et al.' (1994") estimate from ROSAT data that 
that the total luminosity of the Virgo cluster in the 0.1 — 
2.4 keV band is Lx 8x lO^^erg s"^ Only 70% of this power 
comes from the hot gas in the halo of M87, while most of the 
rest comes from M49 and M86. Note that this observational 
estimate is compensated by the fact that another 50% of Lx 
comes from outside the observation band. The final luminos- 
ity of our simulated cluster, Lx ^ 9 x lO^^ergs"^, is there- 
fore very close to value inferred from observations of M87. 
To our opinion, it is even more significant that the black hole 
accretion r ate we have computed is consistent with the value 
inferred bv lPi Matteo et al.1 (|2003l l from studying M87 with 
the CHANDRA X-ray images: they found M. - O.IM© yr"^ 
and LBondi ^ 5 x 10^^ erg s~^. There are also significant mor- 
phological similarities between the cavities opened by the 
jets in o ur simulation and the radio structures observed in 
M87 by (jOwen et al.ll2000l l (see Figure 6). The duration of 
the active phase, tufe ~ 2.5 x 10^ yr, is also in the range of 
what is expected for the typical lifetime of a radio source. 
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APPENDIX A: THE KOMATSU & SELJAK 
HYDROSTATIC SOLUTION FOR THE HOT 
GAS OF AN NEW HALO 

If the distribution of the dark m atter is governed by the 
NFW (|Navarro et al.lll996l . [1997 ) density profile, then the 
equation of hydrostatic equilibrium reads: 



1 dp _ GMvir ln(l + cx) - cx/{l + cx) 
p dr ln(l + c) — c/(l + c) 



(Al) 



Here Mvir is the halo mass within the virial radius r^ir and 
X = r /ryir. For th e NFW concentration parameter, c, we 
follow ISeliakI (|2000l ) and adopt: 

c = 6(Mvir/lO'^M0)-°-^ (A2) 

iBullock et al ]|2001) have argued that a model in which the 
concentration depends on the redshift of collapse agrees bet- 
ter with the observation, but Eq. ()A2|) is perfectly adequate 
for the purpose of constructing a density profile that approx- 
imates that of the ICM of a galaxy cluster. 

To be able to solve Eq. ()A1|) we must introduce an ad- 
ditional constraint on the temperature profile. An isother- 
mal sphere is inconsistent with the requirement that the 
gas trac es the dark matter in the outer parts of the halo. 
Instead, iKoma tsu & Seljak (2001) find a good agreement 
with the measured X-ray profiles and with the observed 
mass-temperature relation when the gas is described as a 
polytrope, with the relation 



P= (7p - l)poeo(p/po)^P. 



(A3) 



Here po and eo are the density and the specific internal en- 
ergy of the gas at the bottom of the potential well, while 
7p is the polytropic index. By substituting Eq. ()A3|) into 
Eq. (ED) , we find that: 



' po y dx po 



1 1 ln(l + cx) - cx/{l + cx) 
'rj^x^ ln(l + c) -c/(l + c) ' 



(A4) 



where 



^0 



GMvir 



Tvir/cTo 



(A5) 

' vir (7p - l)eo rvir/cTo 

and To is the gas temperature in the cluster core. The general 
solution of Eq. ()A4|) is 



po 



7p-1 



7p7yo ln(l + c) - 



1 



ln(l + cx) 



(A6) 



l+c L 

where the parameters po and 770 are determined by imposing 
that p satisfies the conditions: 
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/>(rvir) = /bPdm(rvir) 

and 



dlnp 



(rvir) 



d Inpdn 



(rvir). 



(A7) 



(A8) 



d Inr d Inr 

Here pdm is the dark matter density distribution in the NFW 
model and /b =0.1 is the baryonic mass fraction. Eq. ()A7|) 
gives: 

/bPdm(rvir) 3Mvir 



Po 



PO 

yo(rvir) 



Pvir 



(A9) 



q. ()A8|) gives: 



7p 



c+1 , . 
3^ + ^^^ 



3(l + c)2ln(l + c) -c/(l + c 
c-ln(l + c) 



■Pvir. 



'ln(l + c) -c/(l + c) 



(AlO) 



The conditions ()A7|) and (|A8p guarantee that the gas density 
profile and the dark matter density profile have the appro- 
priate normalization and have the same slope at r rvir, 
but do not guaran tee that the two pro f ile wi ll not diverge at 
r > Tvir- However. iKomatsu SeliakI (|200li) point out that 
one can fine tune the polytropic index 7p to 



7p = 1.15 + 0.01(c- 6.5) 

so that p fhpdm for rvir/2 < r < 2rv 



(All) 



